unit tools;
interface
uses common;
procedure check_lat_center;
procedure mxle3(A:v33;B:v3;Var v:v3);
procedure ortho(x : v2;var y:v2);
procedure ortho3(x : v3;var y:v3);
procedure tra(x : v2;Var y:v2);
procedure tra_perp(x : v2;Var y:v2);
procedure shuffle(mu,ro,mi:integer; Var x:i10);
procedure i5v2(x	: i5;var y:v2);
procedure i5v2p(x : i5;var y:v2);
procedure v5v2(x : v5;var y:v2);
procedure bound5(Var x : i5; Var y:v3);
procedure bound5r(Var x : v5; Var y:v3);
procedure v5v2p(x : v5;var y:v2);
procedure get_atom_basis(ip: atomrec; Var rb:v33);
procedure basis;
procedure add2para(rec : atomrec; ch:integer;w:v3);
procedure add_chem_para(rec : atomrec; ch:integer);
procedure para2xyz(rec : atomrec;ch:integer; ipos,iposd:integer;Var x:v3);
procedure och_status;
procedure para_clean; {remove all chem-parameter entries with zero counts}
function strim(x : mystring):mystring;
function randu(Var ir:longint):real;
function sp(x,y:v2):myreal;
function po(i  : integer):integer;
function dperp_offset(offs :integer; x : integer):integer;
function nui(Var x : i5 ):integer;
function bound(z: myreal):myreal;
function absv(x	: v2):myreal;
function absv3(x : v3):myreal;
function dist(u,v : v2):myreal;
function dist3(u,v : v3):myreal;
function teq10(x,y : i10; var a:i10):boolean;
function eq10(x,y : i10):boolean;
function eq5(x,y : i5):boolean;
function eqatomrec(x,y : atomrec):boolean;
function eqpath(p :integer; x,y : vpath):boolean;
function eqpath012(p :integer; x,y : vpath):boolean;
function teq_node(po,x : integer): boolean;

implementation

procedure och_status;
{sets GLOBALS nnat, naoch, pele, ele, nnele}
Var
   i,j,k,nat,o,p : integer;
   orb,ch	 : integer;
begin
   nat:=0;
   for i:=1 to t_orb do for j:=0 to z_ele do begin
      naoch[i,j]:=0;
   end;
   for j:=0 to z_ele do nnat[j]:=0;
   maxorb:=0;
   for i:=1 to tnodes do begin
      for j:=1 to nan[i] do begin
	 orb:=n_a[i,j,2]; ch:=n_a[i,j,1];
	 for o:=1 to 3 do begin
	    p:=n_para[i,j,o];
	    o3pxyz[orb,o]:=p;
	 end;
	 if (orb>maxorb) then maxorb:=orb;
	 if (ch>0) then begin
	    nat:=nat+1;
	    naoch[orb,ch]:=naoch[orb,ch]+1;
	 end;
	 nnat[ch]:=nnat[ch]+1;
      end;
   end;
   k:=0; for i:=1 to z_ele do if (nnat[i]>0) then begin
      k:=k+1; ele[k]:=i; pele[i]:=k;
   end;
   nnele:=k;
   natoms:=nat;
end; { och_status }
procedure para_clean; {remove all chem-parameter entries with zero counts}
Var
   i,j,k,orb : integer;
begin
   OCH_STATUS;
   for i:=1 to n_pxyz do begin
      orb:=pxyz2orb[i];
      for j:=1 to nn_pxyz[i] do begin
	 if (naoch[orb,ch_pxyz[i,j]]=0) then begin {shift}
	    for k:=j+1 to nn_pxyz[i] do begin
	       ch_pxyz[i,k-1]:=ch_pxyz[i,k];
	       oche[orb,k-1]:=oche[orb,k];
	       dpxyz[i,k-1]:=dpxyz[i,k];
	       nn_pxyz[i]:=nn_pxyz[i]-1;
	       noche[i]:=noche[i]-1;
	    end;
	 end;
	 if (naoch[orb,ch_pxyz[i,1]]=0) then begin
	    nn_pxyz[i]:=0; noche[orb]:=0;
	 end;
      end;
   end;
end; { para_clean }
function strim(x : mystring):mystring;
{get rid of spaces at the word beg: there is pascal }
Var
   j,k : integer;
   y   : mystring;
begin
   y:='';
   j:=length(x);
   k:=0; repeat k:=k+1;
      if ((x[k]<>' ') and (x[k]<>'#')) then y:=y+x[k];
   until (k=j) or (x[k]='#');
   strim:=y;
end;
function randuold(Var ir:longint):real;
var
   xx : real;
   i0  : longint;
begin
   i0:=ir;
   xx:=1048576.0;
   ir := (ir*1025+101) mod 1048576;
   {if ir<=0 then writeln('RND: ',i0:1,' ',ir:1,' ',ir/xx:1:8);}
   {we dont want zero's}
   if ir=0 then i0:=1 else i0:=ir;
   randuold := i0/xx;
end; { randu }
function randu(Var ir:longint):real;
begin
   randu:=random(ir)/ir;
end; { randu }
procedure check_lat_center;
var
   o,i : integer;
   x5  : v5;
begin
   {is the lattice centered?}
   lat_centered:=true;
   for o:=1 to 5 do x5[o]:=0.0;
   for i:=1 to 2 do 
      for o:=1 to 5 do x5[o]:=x5[o]+0.5*bai[i,o];
   {write('LATTICE CENTERING: ');
   for o:=1 to 5 do write(round(2*x5[o]):1,' '); writeln;}
   for o:=1 to 5 do if odd(round(2*x5[o])) then lat_centered:=false;
end; { check_lat_center }
function sp(x,y:v2):myreal;
begin
   sp:=x[1]*y[1]+x[2]*y[2];
end; { sp }
function po(i  : integer):integer;
var j:integer;
begin
   j:=i mod 10;
   if (j<=0) then j:=j+10;
   po:=j;
end;
function dperp_offset(offs :integer; x : integer):integer;
begin
   x:=x+offs;
   repeat
      if (x>2) then x:=x-5;
      if (x<-2) then x:=x+5;
   until abs(x)<3;
   dperp_offset:=x;
end; { dperp_offset }
function nui(Var x : i5 ):integer;
var q,o,y : integer;
begin
   if (dpmodel=0) then nui:=0;
   if (dpmodel=1) then begin {STANDARD: 5 levels in discrete perp space}
      q:=0; for o:=1 to 5 do q:=q+x[o];
      y:=0;
      repeat
	 if (q<-2) then begin q:=q+5; y:=y+1; end;
	 if (q>2) then begin q:=q-5; y:=y-1; end;
      until (abs(q)<3);
      for o:=1 to 5 do x[o]:=x[o]+y;
      nui:=q;
   end;
   if (dpmodel=2) then begin {discrete perp is "even-odd": 01010101....}
      q:=0; for o:=1 to 5 do q:=q+x[o];
      {nui:=0;}
      nui:=-1;
      if odd(q) then nui:=1;
   end;
end; { nui }
procedure mxle3(A:v33;B:v3;Var v:v3);
Var o,i,j:integer;
    det,dx1,dx2,dx3:real;
begin
   v[1]:=0; v[2]:=0; v[3]:=0;
   DET:=A[1,1]*A[2,2]*A[3,3]+A[1,2]*A[2,3]*A[3,1]
   +A[2,1]*A[3,2]*A[1,3]-A[1,3]*A[2,2]*A[3,1]
   -A[1,1]*A[2,3]*A[3,2]-A[3,3]*A[1,2]*A[2,1];
   if (abs(det)<1e-14) then begin
      writeln('MXLE3 error: det=0'); halt;
   end;
   DX1:=B[1]*A[2,2]*A[3,3]+A[1,2]*A[2,3]*B[3]
   +B[2]*A[3,2]*A[1,3]-A[1,3]*A[2,2]*B[3]
   -B[1]*A[2,3]*A[3,2]-A[3,3]*A[1,2]*B[2];
   DX2:=A[1,1]*B[2]*A[3,3]+B[1]*A[2,3]*A[3,1]
   +A[2,1]*B[3]*A[1,3]-A[1,3]*B[2]*A[3,1]
   -A[1,1]*A[2,3]*B[3]-A[3,3]*B[1]*A[2,1];
   DX3:=A[1,1]*A[2,2]*B[3]+A[1,2]*B[2]*A[3,1]
   +A[2,1]*A[3,2]*B[1]-B[1]*A[2,2]*A[3,1]
   -A[1,1]*B[2]*A[3,2]-B[3]*A[1,2]*A[2,1];
   v[1]:=dx1/det; v[2]:=dx2/det; v[3]:=dx3/det;
end; { mxle3 }
function bound(z: myreal):myreal;
begin
   bound:=z;
   if z>0.499999 then bound:=z-1.0;
   if z<-0.500001 then bound:=z+1.0;
end; { bound }
function absv(x	: v2):myreal;
begin
   absv:=sqrt(x[1]*x[1]+x[2]*x[2]);
end; { absv }
function absv3(x : v3):myreal;
begin
   absv3:=sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]);
end; { absv3 }
procedure ortho(x : v2;var y:v2);
var o,i	: integer;
begin
   for o:=1 to 2 do y[o]:=0;
   for i:=1 to 2 do begin
      for o:=1 to 2 do y[i]:=y[i]+x[o]*ba[o,i];
   end;
end; { ortho }
procedure ortho3(x : v3;var y:v3);
var o,i	: integer;
begin
   for o:=1 to 3 do y[o]:=0.0;
   for i:=1 to 3 do begin
      for o:=1 to 3 do y[i]:=y[i]+x[o]*bl3[o,i];
   end;
end; { ortho3 }
function dist(u,v : v2):myreal;
var z,w	: v2;
   o	: integer;
begin
   for o:=1 to 2 do z[o]:=bound(u[o]-v[o]);
   ortho(z,w);
   dist:=absv(w);
end;
function dist3(u,v : v3):myreal;
var z,w	: v3;
   o	: integer;
begin
   for o:=1 to 3 do z[o]:=bound(u[o]-v[o]);
   ortho3(z,w);
   dist3:=absv3(w);
end;
procedure tra(x : v2;Var y:v2);
begin
   y[1]:=(x[1]*ba[2,2]-x[2]*ba[2,1])/(ba[1,1]*ba[2,2]-ba[1,2]*ba[2,1]);
   y[2]:=(x[2]*ba[1,1]-x[1]*ba[1,2])/(ba[1,1]*ba[2,2]-ba[1,2]*ba[2,1]);
end; { tra }
procedure tra_perp(x : v2;Var y:v2);
begin
   y[1]:=sqrt(5)*2*(x[1]*bap[2,2]-x[2]*bap[2,1])/(bap[1,1]*bap[2,2]-bap[1,2]*bap[2,1]);
   y[2]:=sqrt(5)*2*(x[2]*bap[1,1]-x[1]*bap[1,2])/(bap[1,1]*bap[2,2]-bap[1,2]*bap[2,1]);
end; { tra }
procedure shuffle(mu,ro,mi:integer; Var x:i10);
{mu:mult; ro: nx5-fold rot; mi:mirror along [1]}
Var i,j,o,k,rox	: integer;
   y		: i10;
   m5		: i5;
begin
   {mi:=0;}{!!}
   for o:=1 to 10 do y[o]:=0;
   m5[1]:=1; m5[2]:=5; m5[3]:=4; m5[4]:=3; m5[5]:=2; {mirror along [1]}
   ro:=abs(ro) mod 5; {ro taking values [0-4]}
   mi:=abs(mi) mod 2; {mi taking values [0-1]}
   {write('shuffle: ',mu:1,' ',ro:1,' ',mi:1,'; ');
   for o:=1 to mu do write(x[o]:1,' '); write(' -> ');}

   {automatically assume ro=0, mi=1}
   {if (mu=2) then begin
   y[1]:=x[2];
   y[2]:=x[1];
   end;}
   {random choice: better for tile reshuffling}
   if ((mu=2) and (mi=1)) then begin y[1]:=x[2]; y[2]:=x[1]; end
   else begin y[1]:=x[1]; y[2]:=x[2]; end;
   if (mu=5) then begin
      for i:=1 to 5 do begin
	 j:=(i+ro-1) mod 5 + 1;
	 if (mi=1) then j:=m5[j];
	 y[i]:=x[j];
      end;
   end;
   if (mu=10) then begin {first 5 are 5-fold rotated, second 5 their mirror}
      for i:=1 to 5 do begin
	 j:=(i+ro-1) mod 5 + 1; {rotate 1st five}
	 if (mi=1) then j:=j+5; {mirror 1st five}
	 y[i]:=x[j];
	 rox:=-ro+5; {second 5: rotate in opposite sense}
	 j:=(i+rox-1) mod 5 + 6; {rotate: 6->10,7->6...,10->9}
	 if (mi=1) then j:=j-5; {mirror 2nd five}
	 y[i+5]:=x[j];
      end;
   end;
   x:=y;
   {for o:=1 to mu do write(x[o]:1,' '); writeln;}
end; { shuffle }
procedure i5v2(x	: i5;var y:v2);
var o,i	: integer;
begin
   for o:=1 to 2 do y[o]:=0.0;
   for o:=1 to 2 do
      for i:=1 to 5 do y[o]:=y[o]+x[i]*e[i,o];
end; { i5v2 }
procedure i5v2p(x : i5;var y:v2);
var o,i	: integer;
begin
   for o:=1 to 2 do y[o]:=0.0;
   for o:=1 to 2 do
      for i:=1 to 5 do y[o]:=y[o]+x[i]*ep[i,o];
end; { i5v2 }
procedure v5v2(x : v5;var y:v2);
var o,i	: integer;
begin
   for o:=1 to 2 do y[o]:=0.0;
   for o:=1 to 2 do
      for i:=1 to 5 do y[o]:=y[o]+x[i]*e[i,o];
end; { i5v2 }
procedure bound5(Var x : i5; Var y:v3);
var u,v	   : v2;
   o,i,j,k : integer;
   s	   : array[1..2] of integer;
begin
   i5v2(x,u);
   tra(u,v);
   for o:=1 to 2 do s[o]:=round(v[o]+0.00001);
   for o:=1 to 2 do begin
      for i:=1 to 5 do x[i]:=x[i]-s[o]*bai[o,i];
   end;
   for o:=1 to 2 do y[o]:=v[o]-s[o]; y[3]:=0.0;
end; { bound5 }
procedure bound5r(Var x : v5; Var y:v3);
var u,v	   : v2;
   o,i,j,k : integer;
   s	   : array[1..2] of integer;
begin
   v5v2(x,u);
   tra(u,v);
   for o:=1 to 2 do s[o]:=round(v[o]+0.00001);
   for o:=1 to 2 do begin
      for i:=1 to 5 do x[i]:=x[i]-s[o]*bai[o,i];
   end;
   for o:=1 to 2 do y[o]:=v[o]-s[o];
   y[3]:=0.0; {should have discrete-perp here?}
end; { bound5 }
procedure v5v2p(x : v5;var y:v2);
var o,i	: integer;
   r	: v3;
begin
   bound5r(x,r);
   for o:=1 to 2 do y[o]:=0.0;
   for o:=1 to 2 do begin
      for i:=1 to 5 do y[o]:=y[o]+x[i]*ep[i,o];
      for i:=1 to 2 do y[o]:=y[o]-r[i]*bap[i,o];
   end;
   if lat_centered then begin
      if abs(r[1]+0.0001)>abs(r[2]) then begin
	 for o:=1 to 2 do y[o]:=y[o]-(bap[2,o])*sqrt(5)/4;
      end else begin
	 for o:=1 to 2 do y[o]:=y[o]+(bap[2,o])*sqrt(5)/4;
      end;
   end;
end; { i5v2 }
procedure get_atom_basis(ip: atomrec; Var rb:v33);
Var m,n,o : integer;
   rbt	      : v33;
begin
   for m:=1 to 3 do
      for n:=1 to 3 do rb[m,n]:=0.0;
   for m:=1 to 2 do
      for o:=1 to 2 do rb[m,o]:=re[ip[3+m],o]*lpar[1];
   rb[3,3]:=ip[6]*lpar[2];
   for m:=1 to 3 do
      for n:=1 to 3 do rbt[m,n]:=rb[n,m];{transpose}
   rb:=rbt; {copy transposed}
end; { get_atom_basis }
function teq10(x,y : i10; var a:i10):boolean;
{x is connectivity table, y is pattern. y[i]=2 is wildcart,
y[i]=1 must be occupied, y[i]=0 must be unoccupied.
x[i]=0 for unoccupied, and atoms# if occupied.
When y[i]=2, x[i]=-atom#.}
var i,k : integer;
   b  : boolean;
begin
   b:=true; for i:=1 to 10 do begin
      k:=y[i];
      if((x[i]>0) and (k=0)) then b:=false;
      if((x[i]=0) and (k=1)) then b:=false;
      {if (((x[i]>0) and (k=0)) or ((x[i]=0) and (k=1))) then b:=false;}
      if((x[i]=0) and (k<0)) then b:=false;
      if((x[i]>0) and (k<0)) then if(tvl[x[i]]<>k) then b:=false;
      {if((x[i]>0) and (k=3) and (tvl[x[i]]<>-ttvl)) then b:=false;}
      a[i]:=x[i];
      if (k=2) then a[i]:=-x[i];
   end;
   teq10:=b;
end; { teq10 }
function eq10(x,y : i10):boolean;
var i : integer;
begin
   eq10:=true; for i:=1 to 10 do if (x[i]<>y[i]) then eq10:=false;
end; { eq10 }
function eq5(x,y : i5):boolean;
var i : integer;
begin
   eq5:=true; for i:=1 to 5 do if (x[i]<>y[i]) then eq5:=false;
end; { eq5 }
function eqatomrec(x,y : atomrec):boolean;
var i : integer;
begin
   eqatomrec:=true;
   for i:=1 to natomrec do if (x[i]<>y[i]) then eqatomrec:=false;
end; { eqatomrec }
function eqpath(p :integer; x,y : vpath):boolean;
var i : integer;
begin
   eqpath:=true; for i:=1 to p do if (x[i]<>y[i]) then eqpath:=false;
end; { eqpath }
function eqpath012(p :integer; x,y : vpath):boolean;
var i : integer;
begin
   eqpath012:=true; for i:=1 to p do if (y[i]<>2) then begin
      if (((x[i]>0) and (y[i]=0)) or ((x[i]=0) and (y[i]=1))) then eqpath012:=false;
   end;
end; { eqpath012 }
function teq_node(po,x : integer): boolean;
{po: node pointer; ds:top/bottom; x=[0..2] record}
begin
   if (x=2) then teq_node:=true else begin
      teq_node:=((po>0) and (x<>0)) or ((po=0) and (x=0));
      if((po>0) and (x<0)) then teq_node:=(tvl[po]=x);
   end;
end;
procedure basis;
var i,j,k,o : integer;
   luc,fib  : array[0..50] of integer;
   inf5	    : array[1..5] of i5;
begin
   pi:=4.0*arctan(1.0);
   for i:=1 to 5 do
   begin
      {parallel space}
      e[i,1]:=cos(2.0*pi*(i-1)/5.0); 
      e[i,2]:=sin(2.0*pi*(i-1)/5.0);
      {perp space}
      ep[i,1]:=cos(2*2.0*pi*(i-1)/5.0); 
      ep[i,2]:=sin(2*2.0*pi*(i-1)/5.0);
   end;
   {tau-inflation matrix:}
   for i:=1 to 5 do for j:=1 to 5 do inf5[i,j]:=0;
   inf5[1,3]:=-1; inf5[1,4]:=-1;
   inf5[2,4]:=-1; inf5[2,5]:=-1;
   inf5[3,5]:=-1; inf5[3,1]:=-1;
   inf5[4,1]:=-1; inf5[4,2]:=-1;
   inf5[5,2]:=-1; inf5[5,3]:=-1;
   {star re[] for use in the table lt[]}
   for i:=1 to 10 do
   begin
      re[i,1]:=cos(2.0*pi*(i-1)/10.0); 
      re[i,2]:=sin(2.0*pi*(i-1)/10.0);
   end;
   {5d representation of the star re[]}
   for i:=1 to 10 do for j:=1 to 5 do bs[i,j]:=0;
   for i:=1 to 5 do begin
      bs[2*i-1,i]:=1;
      j:=(i+3) mod 5;
      if (j=0) then j:=5;
      bs[2*i,j]:=-1;
   end;
   for i:=1 to 10 do for j:=1 to 5 do infl[i,j]:=0;
   for i:=1 to 10 do begin
      for j:=1 to 5 do for k:=1 to 5 do infl[i,k]:=infl[i,k]+bs[i,j]*inf5[j,k];
   end;
   for i:=1 to 10 do begin
      {for o:=1 to 5 do write(infl[i,o]:3); writeln;}
   end;
   tau:=(1.0+sqrt(5.0))/2.0;
   eta:=sqrt(tau+2.0)/tau;
end; { basis }
procedure add2para(rec : atomrec; ch:integer;w:v3);
{NOTICE: using globals nn_pxyz, ch_pxyz, dpxyz}
Var
   i,j,o,po : integer;
   done	    : boolean;
   mu	    : integer;
begin
   for i:=1 to 3 do if (rec[i]>0) then begin {loop 3 parameters bound to orbit}
      po:=rec[i];
      done:=false;
      mu:=1; if(i=3)then mu:=multi;
      for j:=1 to nn_pxyz[po] do if (ch_pxyz[po,j]=ch) then begin
	 done:=true;
	 {writeln('add2para ',i:1,' ',j:1,' ',po:1,' ',pxyz[po]:1:5,' ',w[i]:1:5);}
	 d_pxyz[po,j]:=d_pxyz[po,j]+mu*w[i]; {mu is z-correction for multi}
	 ndps[po,j]:=ndps[po,j]+1;
      end;
      if not(done) then begin
	 writeln('Add2para: orbit ',pxyz2orb[po]:1,' missing chemistry ',ch:1);
	 halt;
      end;
   end;
end; { add2para }
procedure add_chem_para(rec : atomrec; ch:integer);
{NOTICE: using globals nn_pxyz, ch_pxyz, dpxyz}
Var
   i,j,o,po,p2o	: integer;
   done		: boolean;
begin
   if (ch>0) then begin
      for i:=1 to 3 do if (rec[i]>0) then begin {loop 3 params bound to orbit}
	 po:=rec[i];
	 p2o:=pxyz2orb[po];
	 if (nn_pxyz[po]=0) then begin
	    nn_pxyz[po]:=1;
	    ch_pxyz[po,1]:=ch;
	    dpxyz[po,1]:=0.0;
	    noche[p2o]:=1;
	    oche[p2o,1]:=ch;
	 end else begin
	    done:=false;
	    for j:=1 to nn_pxyz[po] do if (ch_pxyz[po,j]=ch) then done:=true;
	    if not(done) then begin
	       nn_pxyz[po]:=nn_pxyz[po]+1;
	       noche[p2o]:=nn_pxyz[po];
	       if (nn_pxyz[po]>t_ele) then begin
		  writeln('add_chem_para: Number of chemical species exceeds t_ele.');
		  for o:=1 to nn_pxyz[po]-1 do write(ch_pxyz[po,o]:1,' ');
		  writeln(' -> trying to add ',ch:1);
		  halt;
	       end;
	       ch_pxyz[po,nn_pxyz[po]]:=ch;
	       oche[p2o,nn_pxyz[po]]:=ch;
	       dpxyz[po,nn_pxyz[po]]:=0.0;
	       if debug[4] then begin
		  write('ADD_CHEM_PARA: n par orb ch ',i:3,' ',po:3,' ',p2o:3,' ',ch:3,' -> ');
		  for o:=1 to nn_pxyz[po] do write(ch_pxyz[po,o]:3); writeln;
	       end;
	    end;
	 end;
      end;
   end;
end; { add_chem_para }
procedure para2xyz(rec : atomrec;ch:integer; ipos,iposd:integer;Var x:v3);
{NOTICE: using globals nn_pxyz, ch_pxyz, dpxyz}
Var i,j,o,po,ma,ma0,pox	: integer;
   done			: boolean;
   xx			: v3;
   x2,y2		: v2;
begin
   for o:=1 to 3 do xx[o]:=0.0;
   if (ch>0) then begin
   {check we have chemical parameter for given atom}
      ma:=0; ma0:=0; for i:=1 to 3 do if (rec[i]>0) then begin
	 ma0:=ma0+1;
	 for j:=1 to nn_pxyz[rec[i]] do
	    if (ch_pxyz[rec[i],j]=ch) then ma:=ma+1;
      end;
      {NOTE: case ma0=0 for completely fixed orbit is treated automatically}
      {add new chemistry}
      if (ma0<>ma) then begin
	 if (ma>0) then begin
	    writeln('Atom-chemistry parameter error: chemistry cannot be different for different cartezian components!');
	    for o:=1 to 6 do write(rec[i]:1,' '); writeln;
	    halt;
	 end;
	 add_chem_para(rec,ch);
      end;
   end;
   {optional object origin (rec[8]>0)}
   if (rec[8]>0) then for o:=1 to 2 do xx[o]:=xx[o]+ipos*obj_ori[rec[8],o];
   for i:=1 to 3 do if (rec[i]>0) then begin
      po:=rec[i]; {parameter pointer}
      pox:=rec[i+3]; {basis pointer}
      if (i<3) then begin {deca-plane}
	 for o:=1 to 2 do xx[o]:=xx[o]+ipos*pxyz[po]*re[pox,o];
	 {chemistry-specific displacement: ideal position for ch=0}
	 for j:=1 to nn_pxyz[po] do if (ch=ch_pxyz[po,j]) then begin
	    for o:=1 to 2 do xx[o]:=xx[o]+iposd*dpxyz[po,j]*re[pox,o];
	 end;
      end else begin {periodic direction; pox is +-1}
	 xx[i]:=bound(ipos*(pxyz[po]*pox+0.5*rec[7]));
	 {chemistry-specific displacement; for ch=0 nn_pxyz=0}
	 for j:=1 to nn_pxyz[po] do if (ch=ch_pxyz[po,j]) then begin
	    xx[i]:=xx[i]+iposd*dpxyz[po,j]*pox;
	 end;
      end;
   end;
   for o:=1 to 2 do x2[o]:=xx[o];
   tra(x2,y2);
   for o:=1 to 2 do x[o]:=y2[o];
   x[3]:=xx[3]/multi;
   {NOTE: X is now in crystallographic format!!! (without R-node component)}
end;

end.


